Hydrodynamic collective modes for cold trapped gases 



O 

> 
O 

(N 



C/3 



I 

o 



(N 
> 

m 
in 

rn 
O 



X 



Igor Boettcher/ Stefan Floerchinger/'^ and Christof Wetterich-'^ 

^ Institut fiir Theoretische Physik, Universitdt Heidelberg, 
Philosophenweg 16, D-69120 Heidelberg, Germany 
^Physics Department, Theory Unit, CERN, CH-1211 Geneve 23, Switzerland 
(Dated; January 18, 2013) 

We suggest that collective oscillation frequencies of cold trapped gases can be used to test pre- 
dictions from quantum many-body physics. Our motivation lies both in rigid experimental tests of 
theoretical calculations and a possible improvement of measurements of particle number, chemical 
potential or temperature. We calculate the effects of interaction, dimensionality and thermal fluc- 
tuations on the collective modes of a dilute Bose gas in the hydrodynamic limit. The underlying 
equation of state is provided by non-perturbative Functional Renormalization Group or by Lee- 
Yang theory. The spectrum of oscillation frequencies could be measured by response techniques. 
Our findings are generalized to bosonic or fermionic quantum gases with an arbitrary equation of 
state in the two-fluid hydrodynamic regime. For any given equation of state P{^i, T) and normal 
fluid density n„{fJ.,T) the collective oscillation frequencies in a d-dimensional isotropic potential are 
found to be the eigenvalues of an ordinary differential operator. We suggest a method of numerical 
solution and discuss the zero-temperature limit. Exact results are provided for harmonic traps and 
certain special forms of the equation of state. We also present a phenomenological treatment of 
dissipation effects and discuss the possibility to excite the different eigenmodes individually. 

PACS numbers: 03.75.Kk, 67.85.-d, 05.30.Jp 



I. INTRODUCTION 

Ultracold quantum gases provide a powerful exper- 
imental tool to investigate many-body phenomena in 
quantum systems Their distinguished role arises 

from mainly two facts. First they constitute clean realiza- 
tions of microscopic Hamiltonians. For example, whereas 
a continuum gas of fermions with contact interaction or 
the (Bose-) Hubbard-model are only approximations to 
solid state systems, they provide realistic descriptions 
for cold trapped gases. Therefore the latter allow for 
a direct verification of theoretical predictions concerning 
these simple microscopic models. The second important 
feature of ultracold gases is the possibility to tune the 
system parameters by application of external fields and 
to enter the regime of strong interactions with the help of 
Feshbach resonances. Again, this is in clear contrast to 
e.g. solids, where the parameters are fixed by the sample. 

An important and promising direction in cold gases re- 
search is the study of strongly correlated systems. The 
latter arise on all scales of energy, ranging from high- 
temperature superconductors to neutron stars, heavy-ion 
collisions and the phase structure of quantum chromo- 
dynamics. These systems require sophisticated theoreti- 
cal methods, because mean-field theory and perturbative 
techniques break down due to the strong correlations. 
However, since there is no small parameter, the accu- 
racy of non-perturbative approaches is far from being 
well understood. The systematic comparison with ex- 
periments could therefore shed light on the reliability of 
first-principle methods like lattice simulations, the renor- 
malization group, 2PI approaches, extensions of density 
functional theory, conformal field theory etc. 

From our above consideration it should be clear that 



ultracold quantum gases provide an ideal basis for such 
a systematic comparison of theory and experiment. We 
are thus lead to an important question: What observ- 
ables for cold gases reveal advancements and shortcom- 
ings of our theoretical treatment of strongly correlated 
systems? The measurement of these quantities will then 
set stringent bounds on theoretical descriptions and rule 
out insufficient methods. 

In recent experiments, the equation of state has proven 
to be such an observable 0; IM- In this paper we inves- 
tigate to what extend collective oscillation frequencies 
of trapped gases, which can be derived from the equa- 
tion of state, contain similar or additional information. 
For a trapped gas collective modes describe small devi- 
ations from static equilibrium. They typically show an 
oscillatory behavior with characteristic frequencies. The 
number density of atoms in a trap obeys for a particular 
oscillation mode 



n{x,t) — no{x) + Sn{x)cos{iut), 



(1) 



where no{x) denotes the equilibrium density profile and 
we assume small amplitudes, -C no- We will show how 
u! can be calculated for a d-dimensional trap in the ideal 
two-fiuid hydrodynamic regime for an arbitrary equation 
of state. 

For a gas in equilibrium the equation of state consists 
in the pressure P as a function of two independent ther- 
modynamic variables, e.g. chemical potential fj, and tem- 
perature T. In situations where local density approxima- 
tion (LDA) is valid, the equation of state T) can 
be extracted from density images of the trapped cloud 
16;]. We will specify the applicability of LDA below. Here 
we mention only that the system has to be considered 
in equilibrium at each point of the trap separately. This 
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might seem a strong statement at first sight, but, for 
sufficiently high densities and strong interactions, exper- 
iments have shown LDA to be valid at least approxi- 
mately. For example, the coefficient of the Lee-Huang- 
Yang-correction to the equation of state of a bosonic gas 
could be measured to good accuracy [5] . Since the equa- 
tion of state is obtained from a calculation of the par- 
tition function, particular numbers like critical tempera- 
ture and exponents, or even the whole equation of state 
P(/i, T) over a certain range of n and T, can be used to 
verify predictions from many-body theory. If a theoret- 
ical method fails in the context of cold atoms, it might 
also be questionable in other situations. 

For sufficiently low temperatures ultracold quantum 
gases can show superfluidity. Then the equation of state 
consists in two functions, the pressure T) and the 
normal fluid density n„(^, T). Again these quantities can 
be computed from first principles. 

The response of a collision-dominated system in ther- 
modynamic equilibrium to a slight perturbation away 
from equilibrium is completely determined by conser- 
vation of mass, energy and momentum. To first order 
this response can be computed from ideal hydrodynam- 
ics and thus only requires knowledge of the equation of 
state - an equilibrium quantity. Since thermodynamics 
and thus LDA is the static limit of ideal hydrodynam- 
ics, the applicability of the latter might also be given for 
sufficiently high densities and strong interactions. The 
question whether the hydrodynamic limit can be reached 
in experiments is under debate. However, this regime 
is very promising from a theoretical perspective because 
details about the experimental setting become less impor- 
tant. Therefore, already a few quantities like the equa- 
tion of state or transport coefficients can be used to test 
theory. One of the purposes of this paper is to provide an 
idea how measurements in this regime could look like in 
future experiments. If it should turn out to be impossi- 
ble to reach the ideal hydrodynamic regime of cold gases, 
our considerations are still valuable as limiting cases of 
more elaborate treatments. 

In the presence of a non-vanishing superfluid density 
the ideal hydrodynamic equations have to be modified be- 
cause the order parameter enters as an additional macro- 
scopic degree of freedom. The proper two-fluid hydrody- 
namic equations have been derived by Landau. We pro- 
pose that collective modes of trapped gases in this regime 
are an interesting observable in the above sense, i.e. they 
can be used to test theoretical methods which calculate 
the equation of state from a given Hamiltonian. The 
measurement of oscillation frequencies can supplement 
direct determinations of the equation of state from den- 
sity images. In fact, we will show in this paper how the 
spectrum of collective modes can be computed from the 
knowledge of the functions P{^J., T) and nn{^Ji, T) by solv- 
ing an eigenvalue problem. Since thermodynamic deriva- 
tives up to second order enter the calculation, insufficient 
calculations of the equation of state will not yield con- 
sistent results for several modes and over wide ranges of 



external parameters like temperature and particle num- 
ber. Of particular interest is the behavior of the equation 
of state at the second order superfluid phase transition. 

Any observable which is well understood as a func- 
tion of external parameters like temperature, particle 
number or chemical potential can, of course, in turn be 
used to measure these parameters. Since the hydrody- 
namic regime is insensitive to microscopic details, col- 
lective modes can yield information about macroscopic 
properties of the system. Several conditions have to be 
fulfilled for such an investigation to be reasonable in the 
first place: The hydrodynamic regime has to be reached 
in experiment, the equation of state has to be known to 
sufficient accuracy and the mapping between the equa- 
tion of state and collective oscillation frequencies has to 
be understood. Contributing to the latter point is the 
subject of this paper. 

Collective oscillation frequencies have been a promis- 
ing observable since the first experiment on trapped cold 
gases. Therefore, a great amount of literature alre ady 
exists on measurement and calculation Il5l - l29l l 

of these modes for weakly interacting Bose gases, the 
BEC-BCS crossover and other systems at both zero and 
nonzero temperature. While earlier work [sol - fssj has 
partly taken into account additional microscopic prop- 
erties, we concentrate here on the hydrodynamic approx- 
imation. 

The main results of this paper are the following. We 
derive the eigenvalue problem which determines the hy- 
drodynamic modes for a given equation of state P(/x, T) 
and Unil^jT) in a d-dimensional trap. A numerical im- 
plementation based on discretization is put forward. To 
be concrete we consider an isotropic harmonic trap in 
our calculations. However, the extension to arbitrary 
isotropic traps is straightforward, an implementation of 
anisotropic external potentials requires more numerical 
effort but is in principle possible. 

For a zero temperature Bose gas, where hydrodynamics 
is expected to be valid because the system is completely 
superfluid, we calculate the correction of the breath- 
ing mode beyond mean-field in three and two dimen- 
sions. The equation of state P{iJ,) in this case is ob- 
tained from non-perturbative Functional Renormaliza- 
tion Group. We find the expected shift of the breathing 
mode as a function of the gas parameter in three di- 
mensions and provide a similar result in two dimensions, 
where the gas parameter has to be defined differently. 
We emphasize however that our method is not restricted 
to the lowest breathing mode, but rather allows to deter- 
mine the full spectrum of oscillation frequencies including 
higher dipole, quadrupole, etc. modes. 

As a generic example at nonzero temperature we con- 
sider a dilute Bose gas described by Lee- Yang-theory. 
We determine the frequencies over a wide range of tem- 
perature, varying from low temperatures to the normal 
regime above the superfiuid transition temperature. Al- 
though a dilute Bose gas is not necessarily expected to 
fulfill the conditions of hydrodynamics, our method can 
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formally be applied and allows to extract generic features 
that could arise in spectra of oscillation frequencies. 

We discuss the zero temperature limit in the case of 
a Bose gas and find a new set of zero temperature fre- 
quencies due to oscillations of the normal atoms which 
cannot be resolved by purely superfluid hydrodynamics 
of the condensate. We explain why it is important to 
measure several frequencies and how this can be achieved 
with response techniques. 

Due to the narrow constraints of thermodynamics the 
equation of state and thus the collective modes of other 
systems are expected to show a similar behavior. Par- 
ticularly interesting is a strongly interacting Fermi gas, 
where, however, the equation of state is not yet known 
sufficiently, especially below the critical temperature. 
Further experimental and theoretical progress in this di- 
rection will allow for the calculation of the collective 
modes of a unitary Fermi gas as a function of temper- 
ature. 

Our paper is organized as follows. In Sec. |n]we present 
our results on the dilute Bose gas. We are interested in 
contact to experiment and focus on easy readability, giv- 
ing only the relevant equations without derivation. The 
reader will find references to the subsequent sections, 
where derivations of the formulas are presented in a gen- 
eral framework. In Sec. IIIII we derive the eigenvalue 
problem for an arbitrary equation of state at zero and 
non- vanishing temperature. This will yield the formulas 
already used in Sec. |IT1 The numerical implementation to 
obtain the collective frequencies, which is used through- 
out this paper, is given in Sec. IIII Dl We then come to 
our conclusion in Sec. IIVI In App. |X] our chosen sys- 
tem of units is explained and App. [B] summarizes the 
theory of Lee and Yang for the dilute Bose gas. In App. 
[Cl spherical harmonics in d < 3 dimensions are intro- 
duced. Exact results for the zero temperature Bose gas 
are derived in App. |D] Details of our phenomenological 
discussion of response techniques are given in App. |E1 



II. OSCILLATIONS OF A BOSE GAS 

In this section we demonstrate our ideas for the exam- 
ple of a Bose gas. Due to the simplicity of the latter we 
expect our findings to be generic for a broad class of ther- 
modynamic quantum systems. We keep our presentation 
very brief at this point and refer the reader to the next 
section where the general theory of collective oscillations 
is provided in detail. 

At low temperatures and low densities interactions in 
a homogeneous Bose gas can be described completely in 
terms of contact interactions. In this regime the pre- 
cise form of the interaction potential is irrelevant, im- 
plying microscopic universality. In three dimensions the 
associated coupling constant g^o has dimension of length 
and is related to the s-wave scattering length a through 
gsD = 'iirh'^a/m for bosons with mass to. For compos- 
ite bosons all formulas remain valid with small modifica- 



tions. For example, if we aim at describing the weakly 
interacting BEC-side of the BEC-BCS crossover for two 
component fermions, a has to be replaced by the dimer- 
dimer scattering length, which is proportional to the 
fermion scattering length [s^] , and to becomes the dimer 
mass, which is two times the fermion mass. This cor- 
respondence for the equation of state has been nicely 
demonstrated in Ref. [5|. In two dimensions the cou- 
pling constant becomes dimensionless. The effects aris- 
ing in this interesting situation will be investigated later 
in this section. 



A. Three-dimensional dilute Bose gas at zero 
temperature 

For a three-dimensional Bose gas the condition na^ ^ 
1, where n is the density, corresponds to a dilute and 
weakly interacting system. In this case, the mean-field 
result for the energy density at zero temperature is given 
by 0, 0, [3^ {h = ks = m = 1, see App. |A]for units) 



e{n,a) = = 2TTan^. 



(2) 



Taking a derivative with respect to n, we get the chemical 
potential ^ = gson which enters the Gibbs-Duhem rela- 
tion dP = nd/i for the pressure. The equation of state 
written in the grand canonical variables then becomes 



P(M,a) 



(3) 



which is of polytropic type P oc //"+^ with a = 1. At 
zero temperature the system is completely superfiuid and 
will be described by superfluid hydrodynamics. 

The oscillation frequencies of this system in a spherical 
parabolic trap Vext = ^oj'^r'^ are found by solving the 
eigenvalue problem 



Ag{^) 



for the differential operator 



(4) 
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Pi^i^{z) \ dz^ ' ' dz^ 

(5) 

acting on a function g{z). Here d = 3, see Sec. IIIII 
The operator A depends on the equation of state through 
P^in) and P^^(/i), where a superscript denotes differen- 
tiation with respect to fi. The former of these quanti- 
ties corresponds to the density P^ = n while their ratio 
pM/pMM = dn/dP — (? is related to the velocity of sound 
c. Since only this ratio appears in the operator A, the 
prefactor of Eq. ([2]) is of no importance and the frequen- 
cies will not depend on the interaction strength on the 
mean field level. Indeed, Stringari p^j found the analytic 
expression i^ni — (2n^ -|- 2nl -|- 3n -I- /)^/^wo, where n and 
I are integers. Obviously, no thermodynamic conclusions 
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can be drawn from tliis formula and it may at best help to 
determine the trapping frequency when interactions and 
thermal effects are known to be very small. Stringari's 
formula is a special case of 



2n 



1/2 



u}a.n,i= [—ia + n + l + d/2-l) + l] ujq, (6) 



which holds for d-dimensional spherically symmetric, 
harmonic traps and arbitrary polytropic index a. As 
we will show in App. |d1 this formula can be obtained by 
applying a simple polynomial ansatz g{z) = X]fe=o ^fc^*^ 
with z = z/i?^ to Eq. ffl for P (x It is already 

known in the literature [Sg . Ist} . 

Eq. ^ reveals that the possible oscillation frequencies 
of a trapped system form a discrete set. We will see in 
the following that this property is true for an arbitrary 
equation of state P{fi,T). This feature does not arise as 
a result of imposed boundary conditions but is related 
to the domain of definition of the operator A in a differ- 
ent way. We will discuss this issue in Sec. IIII Dl where 
we solve the eigenvalue problem Q on a finite grid by 
discretization. Here, we only remark that the indices n 
and I of Eq. ^ have the same meaning as in the proba- 
bly more familiar case of a quantum mechanical particle 
in a spherical symmetric potential. Similar to the hy- 
drogen atom, n is related to the number of nodes of the 
collective mode, where g{z) plays a role analogous to the 
wave function. Rotation symmetry implies a conserved 
angular momentum I which enters the operator A as a 
parameter. For / = the collective motion is isotropic, 
while for ^ = 1 it is dipolar, and so on. Details can be 
found in Sec. IHII where we derive Eq. ([S]). We expect for- 
mula ([6]) to be applicable only for sufficiently small values 
of n and I, because many nodes or a complicated angular 
structure may enter in conflict with the assumptions of 
hydrodynamics. 

Apparently, for each I > there is a mode with n — 
and frequency ^/Iloq which is independent of the equation 
of state. In particular, the lowest dipole mode is exactly 
at the trapping frequency. Up to now many measure- 
ments focused on the breathing and lowest quadrupole 
mode, n — 1,1 ~ and n — 1,1 ^ 2, respectively. We 
emphasize that the measurement of two, three or more 
frequencies could achieve high precision when combined 
with theory: For a given equation of state an arbitrary 
number of eigenfrequencies can be obtained by our nu- 
merical procedure. Comparison to measurements in well- 
understood parameter regions tests the ability of a given 
theoretical method to compute the equation of state ac- 
curately. On the one hand, it will at least give us an esti- 
mate on the theoretical errors. On the other hand, if the 
agreement is very good, we can take advantage of a non- 
trivial dependence of the frequencies on thermodynamic 
quantities like temperature, density, chemical potential 
or equation of state parameters like a. For example, if we 
predict the lowest lying three monopole modes to have a 
specific temperature dependence, we may conclude from 



measuring these three modes in what temperature region 
we are. 

Motivated by this we are looking for non-trivial re- 
lations between collective modes and thermodynamic 
quantities. We have already seen that a system purely 
described by mean field theory does not show such a be- 
havior. We expect the mean field picture to be valid for 
very small gas parameter na'^ <C 1. However, as this 
parameter increases, higher order interaction effects be- 
come relevant in the equation of state. The leading order 
correction to the ground state energy density ^ at zero 
temperature has been calculated by Lee, Huang and Yang 
and is found to be i38|] 



£LHY{n,a) = 2Tian 1 
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[no:') 



3U/2 



(7) 



The corresponding pressure as a function of chemical po- 
tential to the same order in perturbation theory is given 
by 



Stto 



157r- 



5/2 



(8) 



Apparently, the interaction strength can no longer be 
eliminated by a rescaling of such that the ratio 

pMy'pMM entering the operator A will always depend on 
a. Thus, the frequencies will depend on the coupling 
constant. While for a homogeneous system the thermo- 
dynamic observables are constant in space, for a confined 
gas one usually refers to their values at the center of the 
trap. Denoting the density in the center of the external 
potential by n(0), a'^n(O) provides a small parameter and 
one expects a shift of the breathing mode {n — 1,1 — 0, 
i.e. UJB = y5tJo) [11113 
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[aMO)] 



1/2 



(9) 



as compared to the Stringari formula above. We con- 
clude that we may use the shift for small values of n{0)a^ 
in order to determine either ri(0) or a very precisely if 
the corresponding second quantity is known from another 
method. 

From this example we can deduce a method to make 
high precision measurements using collective modes. If 
we assume the shifts of the frequencies uin,i to be con- 
tinuous in n{0)a^ or /x(0)a^, where /x(0) is the chemi- 
cal potential in the center of the trap, then for small 
gas parameters, n{0)a^ <C 1 or /i(0)a^ ^ 1, the shift 
will be proportional to some power of the gas parame- 
ter. In the case of Eq. © we found this power to be 
1/2. However, this will depend on the system under con- 
sideration. Nevertheless, driving the system through a 
certain range of the gas parameter, e.g. by the use of 
a Feshbach Resonance, and then collecting the result in 
a double-logarithmic plot one will find the results to lie 
on a straight line for small gas parameter. This simple 
scaling behavior can be used for very precise measure- 
ments of n{0)a^ after a proper calibration. The accuracy 
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FIG. 1: Shift of the mean-field breathing mode cjs = \/5i^o for 
the zero temperature Bose gas in three dimensions. We show 
predictions for different methods: Lee-Huang- Yang-formula 
(|8]), Wu-correction (|f 01) with B — 8 and equation of state from 
Functional Renormalization Group 41] (FRG) . The solid line 
corresponds to Eq. (|9]), a is the scattering length and n(0) 
the density in the center of the cloud. 



of this method is directly related to the accuracy in mea- 
suring the frequency shifts. In Fig. [l]we show the shifts 
beyond mean field due to the LHY-correction which we 
obtained by our numerical implementation described in 
Sec. IIIIDI We reproduce the prediction to a high pre- 
cision in the regime of small gas parameter, where it is 
valid. Although this plot only shows the breathing mode 
we emphasize that it is possible to calculate the shifts of 
the whole frequency spectrum which are found to show 
a similar scaling behavior with ^y n{0)a^. Fig. [T] shows 
that an experimental accuracy for frequency shifts bet- 
ter than 10~^ would be useful for exploring the regime 
of small gas parameters. 

The calculation of higher order corrections to the LHY 
equation of state ([7]) is a difficult task. The next-to- 
next-to-leading order has been derived by Wu [i^, 
for hard-sphere bosons. The energy density receives two 
additional terms 



It does not rely on the expansion in a small parameter 
and therefore the most striking effects are expected in 
strongly interacting systems. For more information on 
the FRG see Refs. [43|. It is still under debate whether 
the regime a — >■ oo can be reached for a Bose gas at very 
low temperatures. One expects the condensate to get 
unstable then because of increasing importance of three- 
body losses. 

In Fig. [T] we show the frequency shift obtained from 
the equation of state calculated with FRG at zero tem- 
perature [4l|. We see that higher values of the gas pa- 
rameter n(0)a^ can no longer be identified with a unique 
6ujB/tOB- This is a fully non-perturbative effect. We re- 
call that one of the assumptions in the beginning was 
that the interaction can be approximated to be point- 
like. However, this is only valid if the scattering length a 
is much larger than the microscopic distance where 
one can resolve the details of the interaction potential. 
For cold bosonic atoms this length is given by the typi- 
cal range of the Van der Waals potential. An effective 
ultraviolet cutoff A is necessary in any field theoretic 
treatment of the system where the contact interaction 
appearing in the Lagrangian of the non-relativistic sys- 
tem is not renormalizable in three dimensions. Therefore, 
the introduction of this scale is no relict of approxima- 
tions but has a clear physical meaning. In perturbation 
theory one assumes A to be infinity and thus the equa- 
tion of state of the homogeneous system only involves 
one dimensionless parameter, the gas parameter na^. A 
proper treatment has to account for the appearance of an 
additional length scale A~^ such that two possible combi- 
nations, na^ and aA, describe the system. As a result the 
shift of the breathing mode should rather be plotted as a 
two-dimensional surface depending on these two param- 
eters. For small gas parameter the additional parameter 
aA does not play a role - if we take equations of state 
P{fi,a) for different values of a and then vary n(0), all 
obtained shifts will lie on a line in the double-logarithmic 
plot. This is the expected scaling behavior. For higher 
densities it makes a difference whether we calculate the 
frequencies from a, Ai) or P{fi, a, A2) with Ai 7^ A2, 
as reflected by the spread of the points in Fig [TJ 



£lhy + 2TTan^ 8 



(log {na^)+B). 



(10) 

Neglecting the logarithmic term, this leads to a correction 
which is proportional to /i'^a in Eq. ([5]). In Fig. [T] we 
have set B — 8, see e.g. Ref. [5] for an overview over the 
expected values of B. Obviously, the general behavior of 
the frequency shifts is not influenced. 

For large values of the gas parameter perturbation 
theory is no longer applicable and more sophisticated 
methods are necessary. The Functional Renormalization 
Group (FRG) for the effective average action 4J] is a 
non-perturbative quantum fleld theoretical method and 
is in particular able to describe the full thermodynamics 
of systems which are realized in cold atom experiments. 



B. Lower-dimensional systems 

With regard to formula ([B]) we may wonder whether the 
cases d — 2 and d ~ 1 refer to highly anisotropic traps. 
These are of great relevance for experiments where one is 
often not dealing with a spherical symmetric potential. 
For example, the early experiments in the JILA group 
were performed in a disk-shaped conflnement ^ while 
the MIT group used a cigar-shaped trap [8] . 

A really lower-dimensional system is obtained in a trap 
where the quantum gas is in its ground state in one or 
two directions. This will be the case for very tight con- 
flnement. When calculating the equation of state one can 
then neglect quantum and thermal fluctuations in these 
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FIG. 2: Shift of the breathing mode relative to ljb = 2cjo for 
a two-dimensional dilute Bose gas for different values of the 
coupling constant g2D- We show the dependence on /i(0)A~^, 
with chemical potential fj.{0) in the center of the trap and 
effective UV cutoff A. The equation of state is provided by 
Functional Renormalization Group. 

directions. The mean field Bose gas in this scenario is 
described by Eq. ([B]) when inserting d = 2 or d = 1 
and a ~ 1. In these cases the system is isotropic in a 
lower-dimensional geometry. 

A different situation arises if the system is not in its 
ground state in the directions of tight confinement and 
thus still three-dimensional or in a crossover between 
three and lower dimensions. Formula ^ cannot be ap- 
plied in this case because the assumption of spherical 
symmetry is not justified. The solution of the hydrody- 
namic equations for these anisotropic cases is more in- 
volved. 

We now consider the two-dimensional Bose gas. Its 
mean field equation of state is fi = g2Dn and the fre- 
quency of the breathing mode is found from Eq. ^ to 
be ujb — 2cjo. The equation of state beyond mean field 
can be obtained by a Functional Renormalization Group 
approach (4l| . The coupling constant (72D, which satisfies 
M = g2Dn for small g2D, is dimensionless in two dimen- 
sions and shows a logarithmic running with the physi- 
cal scale on which the experiment is performed. It will 
vanish for an infinitely large system. However, realistic 
experiments are always performed in traps so that in a 
harmonic potential the oscillator length £0 = yh/mw^ 
provides the largest possible length scale of the physics 
under considerations. Therefore, when calculating the 
equation of state one only has to include quantum fiuc- 
tuations with momenta bigger than £q ^, which acts as 
an infrared cutoff. This is also the reason why Bose- 
Einstein condensation is observed experimentally in two 
dimensions for small temperatures Q < T < Tc- For an 
infinitely extended system the long-range order would be 
destroyed by fluctuations for all non- vanishing tempera- 
tures as required by the Mermin- Wagner theorem [46]. 

In Figs. and |3] we plot the shift of the breathing 



FIG. 3: The same setting as in Fig. [2] but for higher values of 
the coupling g2D- We find rather large negative shifts, which 
have a minimum. 

mode corresponding to lob = Swq for the two-dimensional 
Bose gas at zero temperature with equation of state 
-P(Mj52d) from Functional Renormalization Group cal- 
culations [4l|. Since the coupling constant is dimension- 
less there is no gas parameter for such a system. As 
we mentioned already in the three-dimensional case one 
always has a physical ultraviolet cutoff A when dealing 
with a contact interacting in d > 2 dimensions. Since A 
has dimension of inverse length, the dimensionless vari- 
able involving the chemical potential (or similar for the 
density) is ^(0)A~^, where /i(0) is the chemical potential 
of the gas in the center of the trap. A good choice for 
A~^ is the oscillator length of tight trapping. We observe 
from Fig. [2]that the frequencies for small interactions de- 
pend only weakly on /i(0)A~^. For larger coupling g2D 
we find deviations from this behavior in Fig. |3l We em- 
phasize that the dependence on the microscopic cutoff 
for two-dimensional bosonic systems is not a shortcom- 
ing of approximations, but a physical effect. It is due 
to the running of couplings. While the relevant coupling 
is dimensionless, the equation of state retains memory 
of scales in form of a dependence on the physical effec- 
tive cutoff. This effect becomes substantial for large g2D- 
Therefore, the strong dependence of the frequency shift 
on the dimensionless combination /i(0)A^^ for large inter- 
action strength g2D ^ 1 in Fig- [3] signals the breakdown 
of universality. More microscopic properties beyond the 
single interaction parameter g2D are needed to describe 
a concrete experimental realization in this regime. 

We arrive at two interesting experimental scenarios to 
be investigated. On the one hand, by measuring the col- 
lective modes one can distinguish whether one is work- 
ing with a system which is still three-dimensional (disk, 
cigar) or a truly lower dimensional system. On the other 
hand, if the latter regime is reached experimentally it is 
of course tempting to verify the predictions of Fig. |31 
The rather large negative frequency shifts can both test 
theoretical predictions like the occurrence of a minimum. 
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and determine /i(0) or the corresponding density n(0). 



C. Three-dimensional dilute Bose gas for 
non-vanishing temperature 

We now extend our considerations to non-vanishing 
temperature and allow for fluctuations of both density 
and temperature. If we stay below the critical tempera- 
ture of Bose-Einstein condensation the system possesses 
both non-vanishing superfluid and normal fluid density 
and its hydrodynamic behavior has to be described by 
two-fluid hydrodynamics [47l - l49j . 

As explained in Sec. IIIH we again find an eigenvalue 
problem which has to be solved in order to get the col- 
lective frequencies of the system. It has the general form 



A B 
C D 



9i^) 
h{z) 



9{z) 
h{z) 



(11) 



where the differential operators A, B, C and D are de- 
fined in Eqs. - ([7D|. These operators depend on 
the equation of state T) and the normal fluid den- 
sity n„(/i,T). Both quantities have to be provided by 
an underlying microscopic theory. The structure of the 
eigenvalue problem is similar to the simple zero tempera- 
ture case. We will show below that the zero temperature 
limit of Eq. (fTTj) is indeed given by Eq. ^ and this 
behavior will also be found in the frequencies. 

We see that calculating the frequencies at non- 
vanishing temperature is as straightforward as it was at 
T = 0. However, results on the temperature dependence 
of hydrodynamic collective modes are rare in the liter- 
ature [2l|. This is because the equation of state and 
normal fluid density are in general not known or only 
in certain ranges of temperature and chemical poten- 
tial. But when applying the local density approximation 
/i — )• /i — Vexti'r) one drives through a wide interval of 
values for the chemical potential and therefore a com- 
plete equation of state in /i has to be provided. Our mo- 
tivation is in applying the Functional Renormalization 
Group, which is capable of providing the full thermody- 
namics of cold quantum gases, especially the functions 
P(^,T) and n„(/i,T). For the three-dimensional weakly 
interacting Bose gas the equation of state has been calcu- 
lated by Lee and Yang [SO] ■ Their formulas are valid for 
small gas parameter na^ for temperatures not too close 
to the critical temperature. We comment on the criti- 
cal behavior later in this section. The Lee- Yang equa- 
tion of state is summarized in App. [Bj We neglect the 
next-to-leading order LHY-correction ([7]) because we are 
interested in thermal effects. 

In [2l| Shenoy and Ho have calculated the tempera- 
ture behavior of collective modes of a trapped Bose gas 
from two-fluid hydrodynamics and Lee- Yang theory for 
0.6 < T/Tc < 1.2. The authors considered the range of 
parameters where in addition to na^ ^ 1 also anX^ <C 1 
holds. (At = {2Trh'^ /rnksTy^'^ is the thermal wave- 
length.) In this limit the integral expressions for the 




T/T„=0.5 ■ 
Tn"„=0.75 - 



FIG. 4: Equilibrium density profile as a function of the radius 
obtained from Lee- Yang theory using local density approxi- 
mation for a spherical symmetric potential. We observe a 
peak in the center of the trap. In the outer regions the gas is 
in its normal phase. 



equation of state from Lee- Yang theory can be evaluated 
analytically. We will, however, keep the full expressions 
for the equation of state. 

Before discussing the collective modes of a Bose gas 
at non-vanishing temperature we comment on some as- 
pects of the static configuration of the trapped gas which 
are necessary for the interpretation of our results. For a 
homogeneous system with temperature T we can calcu- 
late the critical chemical potential fJ-c{T) of the superfluid 
phase transition. Within the local density approximation 
in the trapped system there will be a radius Tc corre- 
sponding to the phase boundary between the superfluid 
and the normal regions of the cloud. This radius fulflUs 
He = m(0) ~ Vextirc) with the chemical potential /z(0) in 
the center of the trap. The characteristic picture of the 
density profile consists in a narrow condensate peak in 
the center of the trap which is surrounded by a broad 
thermal cloud of the normal gas. Of course, there is also 
a non-vanishing contribution of the normal component 
to the inner regions. We visualize this behavior in FigU 

It is now apparent that for the description of a trapped 
gas at any nonzero temperature one has to know the 
equation of state for both the superfluid and normal 
phase. In particular the presence and dynamics of the 
normal gas are of importance. There can be oscillations 
of the thermal cloud itself and, furthermore, it provides 
a non-trivial background for the oscillations of the con- 
densate. As we approach zero temperature the thermal 
cloud vanishes. Condensate oscillations correspond to so- 
lutions 6n of the hydrodynamic equations which are only 
nonzero inside a sphere with radius Tc- 

The critical temperature Tc of the trapped gas is de- 
fined as the temperature where the condensate peak ap- 
pears in the center of the cloud. This can be reformulated 
as fj,(0) = He or equivalently Tc = 0. For Lee- Yang theory 
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FIG. 5: Oscillations frequencies of collective modes {I = 0) of both the condensate and thermal cloud in an isotropic harmonic 
trap. The equation of state is provided by Lee- Yang theory. In this plot the particle number is fixed to be A'^ = lO'' for each 
value of T/Tc- The oscillations of chemical potential and temperature which correspond to the black solid dots are shown in 
Figs. [B]and[3 



the critical temperature is given by 



= 27r 



V2.9C(3/2) 



(12) 



If we set /i(0) = /ic and calculate the number of par- 
ticles N in the trap, we find the ideal gas prediction 
T^°\n) = fia;o(iV/C(3))i/3 to be satisfied for small gas 
parameter. There are several possibilities to change the 
parameter T/Tc- For fixed /i(0), the particle number will 
change when decreasing the temperature. Another way 
of going to low temperatures is to keep the particle num- 
ber fixed. Then the chemical potential in the center of the 
trap has to be adjusted appropriately. If non-destructive 
measurements of frequencies become possible, a partic- 
ular promising experimental way would be a change of 
T/Tc by a smooth variation of the trap frequency, keep- 
ing the number of trapped atoms fixed. 

In Fig. [5] the isotropic {I — 0) temperature de- 
pendent oscillation frequencies in an isotropic harmonic 
trap are shown for fixed particle number = 10^ and 
a/£o = 0.0005, where £o = ^hjmujQ is the oscillator 
length. This particular choice of parameters implies a 
small gas parameter na'^, which is required in order to 
use Lee- Yang theory, but does not guarantee the ap- 
plicability of hydrodynamics. Nevertheless, we expect 
this approach to capture the characteristics of the tem- 
perature dependence of hydrodynamic collective modes. 
Due to thermodynamic constraints the equation of state 
for the strongly coupled case will share common features 
with Lee- Yang theory and thus will show a qualitatively 
similar oscillation spectrum. For a quantitatively precise 
description of the collective modes of the weakly coupled 
system, which is not the purpose of this paper, one needs 
to include corrections to the hydrodynamic behavior by 
means of coUisionless dynamics. Our results then consti- 
tute a limiting case. 

We observe a rich spectrum of frequencies. The os- 



cillations can be classified as condensate oscillations and 
oscillations of the thermal cloud. The condensate oscil- 
lations correspond to the branches which disappear at 
Tc- The computability of frequencies above the critical 
temperature by our method is a manifestation of the fact 
that the Landau two-fluid model remains formally valid 
for vanishing superfluid component. Of course, the fre- 
quencies for T >Tc could be computed equivalently with 
a one-fluid model for the thermal liquid. We find the os- 
cillation frequencies of the thermal cloud to be very close 
to the spectrum which one obtains using the equation 
of state for an ideal Bose gas and vanishing superfluid 
density. In particular, the lowest mode which also exists 
above Tc is given hy uj — 2loq, see for example Ref. [18]. 

For low temperatures {T/Tc ^ 0.1) the two-fluid re- 
gion, defined hy r < Tc with rc from above, starts to 
grow considerably and the normal fluid in the outer re- 
gion vanishes. However, we were not able to obtain the 
eigenfrequencies of this crossover numerically in a suffi- 
ciently accurate manner. In order to get an understand- 
ing of the very low temperature limit of the collective 
modes we therefore used an implementation where the 
cloud is cut off at rc and thus only the two-fluid region is 
left. One then obtains a spectrum of the type shown in 
Fig. ISl which shows a smooth limit for T — 0, as is dis- 
cussed below. (For increasing T there remains a system- 
atic error in Fig. [8] because with increasing temperature 
the normal mantle surrounding the two-fluid core of the 
cloud cannot be neglected anymore.) 

From Fig. [S] it is clear that there are level-crossing 
frequencies over wide ranges oiT/Tc. These branches can 
be used as a thermometer by measuring a certain number 
of frequencies and then locating them in the plot. How 
this can be done experimentally by response techniques 
is explained in Sec. Ill Dl 

An apparent feature of Fig. [5] is the appearance of 
bifurcations points where one branch of collective exci- 
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tation frequencies splits up into two distinct branches. 
The meaning of these points is very intuitive. While 
for low temperatures superfluid and normal fluid com- 
ponent are closely connected, the modes split up at in- 
termediate values of T/Tc- One of them represents the 
superfluid oscillation, which vanishes at the critical tem- 
perature, the other one corresponds to the oscillation of 
the normal fluid which remains present also above Tc- 
The black solid dots in Fig. [5] correspond to points be- 
fore and after a splitting. We indeed find that for low 
temperatures there is only one solution J/i and 5T for 
this particular eigenvalue, which then smoothly splits up 
into a set of two solutions with different properties. We 
make this statement clear in Figs. H] and [7] where the 
(isotropic) oscillations of chemical potential and temper- 
ature are shown for the frequencies corresponding to the 
black dots in Fig[5l distinguished by the labels "upper" 
and "lower" branch, meaning higher and lower frequency, 
respectively. 

One has to keep in mind, that the oscillations of parti- 
cle number density and entropy density, 5n and 5s, vanish 
at infinity because the thermodynamic functions in Eqs. 
(1551) and ([57| vanish. However, the oscillations of density 
and entropy density are related to 5/i and 5T via 

where the matrix contains the higher derivatives of the 
equation of state P(/x, T) applied to local density approxi- 
mation. (We use here the notation of Sec. IIIIH Since this 
matrix vanishes at infinity, 5^ and 5T can be nonzero. At 
zero temperature, where we have a sharp cloud radius of 
the condensate, the fluctuations of the chemical potential 
(5/i are discontinuously going to zero at this radius. 

From our discussion of the equilibrium density profile 
we see that for all temperatures below the critical temper- 
ature Tc there is a radius Tc where the critical equation of 
state P{fJ-c, T) has to be known. The vicinity of the point 
He might be small but the behavior of the thermody- 
namic functions in this interval could nevertheless influ- 
ence the frequencies significantly. Indeed, the superfluid 
phase transition is of second order and we expect the spe- 
cific heat at constant pressure Cp = T{dS/dT)p^]\[ and 
the isothermal compressibility kt = —V^^{dV/dP)T,N 
to be singular at /ic. Their behavior in the critical re- 
gion as a function of temperature is Cp ^ \T — Tc|~" 
and Kt ~ \T — Tc\~'' with the critical indices 'sH] of 
the three-dimensional XY-universality class. The spe- 
cific heat shows a cusp while the isothermal compress- 
ibility diverges. A systematic study of these effects on 
the collective modes would be very interesting. Lee- Yang 
theory does not allow for such an investigation. 

Using the Lee- Yang formulas given in App. |B]we can 
calculate numerically the coefficient functions appearing 
in the operators A, B, C, and D in Eq. pT|) . For very 
low temperatures we observe that B and C vanish while 
A approaches its zero temperature expression, given by 



in^ =0.215 

T/Tj,=0.25, upper branch 

T/Tj,=0.35, upper branch 
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FIG. 6: Oscillation profile of the chemical potential for the 
frequencies corresponding to the the black solid dots in Fig. 
[S] For T/Tc — 0.215 there is only one solution, which then 
splits up into a lower branch indicated by the dashed line, 
and an upper branch represented by the dot-dashed line. The 
thinner lines correspond to higher temperatures. We observe 
the lower branch to form a bump in the region of the conden- 
sate, while the upper branches smoothes out this boundary. 
The eigenfunctions are normalized here to be 1 in the center 
of the trap. 
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FIG. 7: The same setting as in Fig. [6] but showing the oscil- 
lation profile of temperature. Starting from a solution which 
possesses a cusp, the lower branch pronounces this peak while 
the oscillations corresponding to the upper branch become 
flat. The eigenfunctions are normalized here such that they 
have value 1 at the peak. 

Eq. dSD, 

f A B \ T^o f A{T = 0) \ , . 

[C D ) ^ { D{T^O) )■ 

This corresponds to a decoupling of superfluid and ther- 
mal degrees of freedom. This is expected to happen in 
the zero temperature limit since the few atoms in the 
thermal cloud have no influence on the dynamics of the 
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FIG. 8: Zero temperature limit of the collective oscillations. 
Here, the chemical potential is kept fixed at /x(0)/fiajo = 10, 
which does not coincide with the choice of parameters in Fig. 
[5l The solid lines over the whole range of temperatures cor- 
respond to formula ^ for superfluid oscillations at T = 0. 
The frequencies given by Eq. ((T7|) are indicated by solid lines 
only for small temperatures, where this formula is expected 
to be valid. In this figure, there is a systematic error for 
high temperatures, because the thermal cloud is not treated 
appropriately. 



condensate. The operator entering the eigenvalue prob- 
lem is now block diagonal and all zero temperature fre- 
quencies are reproduced. They correspond to Stringari's 
mean field formula given above. (Note that we have ne- 
glect the LHY-correction.) It is interesting that D does 
not vanish. Indeed, the coefficient functions oi, di and 
^2 entering A and D from Eqs. (|68|) and (|70|) satisfy 



fll T-^0 1 , T^O 1 

^ 3' ~6- 



(15) 



This can be seen by evaluating the integrals for P(/x, T) 
and n„(fj,,T) numerically or applying an expansion for 
T <^ fj, using only the phonon part of the spectrum, see 
App. |D] It can be shown that ai/di corresponds to the 
ratio of the squared first and second velocities of sound, 
and C2, respectively. Therefore we find the well-known 
relation 
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(16) 



to be satisfied for very low temperatures. 

We show in App. |D] that the behavior of the coeffi- 
cient functions of D as T — 0, Eq. ([T5|) . can be used to 
derive the eigenvalues of D by slight modifications from 
the eigenvalues of A, which are given by Eq. ^ with 
a = 1. The corresponding eigenfrequencies are found to 
be 



4n{n + 1) - I 
6 



1/2 



t^o, ("- > !)■ 



(17) 



The experimental verification of this formula will be dif- 
ficult, however. Although for low T (and T — >■ 0) these 



frequencies are solutions to the two-fluid hydrodynamic 
equations, there are only very few atoms left in the nor- 
mal phase to oscillate. From Eq. (fT7|) it is apparent that 
for n = 1, Z = one obtains ^j2/'iujQ which is below the 
trapping frequency. This cannot be achieved by formula 
© with an effective polytropic index a^ff- The exis- 
tence of such modes below the trapping frequency might 
be a special feature of superfiuid hydrodynamics just as 
the appearance of a second velocity of sound. Indeed, in 
their early experiments Stamper-Kurn et al. [8] used a 
cigar-shaped harmonic potential with trapping frequency 
woz/27r — 18.04(1) Hz in axial direction. They observed 
an out-of-phase oscillation between condensate and ther- 
mal cloud of a Bose gas at T = IjiK in the hydrodynamic 
regime. The critical temperature was around 1.7fiK. 
The measured frequency of axial motion was found to 
be ojo/Stt = 17.26(9) Hz, which is below luqz. 

In Fig. [5] we show our numerical results for the fre- 
quencies in the zero temperature limit using the full 
Lee- Yang equation of state for fi{0)/hajo — 10 fixed and 
a/£o = 0.0005. We find formulas © and ^ to be in 
perfect agreement with our data. 



D. Measuring the spectrum with response 
techniques 

We have seen that information is contained both in 
the precise values of single modes and in the whole spec- 
trum of the oscillating system. While the lowest lying 
modes can be excited by varying the parameters of the 
trapping potential (e.g. the potential depth) and then 
measuring the frequency in a free expansion, the whole 
spectrum could be obtained by an other method. In prin- 
ciple one can fit a superposition of arbitrary many modes 
to the freely expanding cloud but the contribution of the 
higher modes will only be weak. We suggest to use a 
response measurement instead to locate the frequencies 
of the trapped cloud. In order to keep the notation short 
we restrict this discussion to the zero temperature case. 

The eigenfrequencies and corresponding eigenfunctions 
obtained so far represent normal coordinates of an os- 
cillating system, the trapped cold gas. The equation 
of motion d^S^x + Ad^j. — has solutions S^nim — 
9niir)r^Yime^^'^"'* which are harmonic oscillations with 
frequency ujni. (More precisely, one should use the oper- 
ator E defined below Eq. (jE2[) instead of A. This is of no 
relevance here, because they have the same eigenvalues.) 

To describe an experimental situation we need to in- 
clude dissipation effects into our description. We there- 
fore introduce a phenomenological damping term by re- 
placing — + rdt, with damping constant F. The 
new eigenfrequencies of the system will be denoted by 
f2„;. We use Fourier decomposition and write Sfi{t) — 
J dr25/i(r2)e~'^*. Since we already know the eigenvalues 
of A we immediately obtain the quadratic equation 



(- 



nl 



(18) 
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for possible solutions The latter simplify due to the 
fact that t2 — 1/r represents a characteristic time scale 
over which damping takes place. In order to observe 
oscillations this scale should be substantially larger than 
ti = l/iOni- The other choice would correspond to the 
overdamped case. Hence lu'^i ^ and we have 



^nl — w„; 1 



iT 

y 



(19) 



for the new eigenfrequencies of the system. Due to damp- 
ing effects they acquire an imaginary part and their real 
part is slightly shifted. This phenomenological treatment 
gives no explanation of the microscopic nature of the 
damping terms. It would be interesting to derive the 
relation between the damping constants F and the five 
transport coefficients rj, ktXij C2 and ^3 of two-fluid hy- 
drodynamics j49| . In principle, one expects F to depend 
on the mode in the form of some continuous function 
F(r2). For simplicity we use here the same phenomeno- 
logical F for all modes. 

The picture of the trapped gas as an oscillating system 
will provide us with an intuition how individual (and thus 
also higher lying) modes can be addressed by virtue of 
response techniques. Consider a classical point particle 
with Hamiltonian H = -^p^ + ^ur^x^ . With a damping 
term included its equation of motion will have the form of 
Eq. (ITS)) and is therefore independent of the mass of the 
particle. However, when applying a small external force 
f(t) the equation of motion for the Fourier components 
of x{t) becomes 

{-^f -iTn + u'^)x{n) = —f{n). (20) 

m 

Obviously, the response x(^) = x{Vt)lf{Vt) is propor- 
tional to the inverse of the mass m. It is therefore more 
difficult to excite a heavy particle than a lighter one. 
We introduce an analog quantity ninim which represents 
the "mass" for the collective oscillation SjjLnim- We call 
iTT-nim the "rcspousc coefficients" . The relation between 
the amplitude of the oscillations of the atom cloud and 
the oscillation of the external perturbation f{ft) is gov- 
erned by rrinim- 

Let us consider a periodic perturbation of the harmonic 
trapping potential, 

Ve^tix) ^ Ve^tix) + cos{nt)FiryYim. (21) 

The perturbation is assumed to be small in comparison 
to Vext- As we show in App. |E]this small variation acts 
as a driving force on the oscillating system psp . There 
it is also shown that in order to excite an oscillation with 
angular dependence (5 oc Yim, SVext has to have the 
same angular dependence. In Eq. (jE12p we give a general 
formula for the response coefficient mnim of the collective 
mode. Given ~ ujni we then have 
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FIG. 9: Response of the zero temperature mean-field Bose gas 
to a perturbation SVext = r''cos{^lt) of the external potential. 



Imx(i^) 



1 



(22) 



for the imaginary part of the response function defined in 
Eq. (|E11[) . The latter is related to dissipation of energy 
in the system. 

For a demonstration we restrict ourselves to the case of 
a three-dimensional mean- field Bose gas, i.e. P(/i) cx: /i^, 
and consider only the isotropic modes with I — 0. In 
addition we assume a monomial perturbation F{r) — r'^ 
with g > 0. This includes the cases of a spatially con- 
stant background field 5Vext{x,t) — cos(ili) {q ~ 0) 
and a time-variation of the trapping frequency, SVext — 
r^cos{ilt) {q — 2). Interestingly, we find in App. lElthat 
not all values of q can be used to excite the whole spec- 
trum. For example, there is no response at all to g = 0, 
while for q — 2 only the breathing mode will be excited 
and for (7 = 4 only the lowest two modes oscillate. We 
show this behavior in Fig. |9l (In our intuitive picture de- 
veloped above the absence of response corresponds to an 
infinite mass.) The normalization in Fig. [S]is such that 
the lowest mode always has the same amplitude. The rel- 
ative amplitudes, both for oscillations with different fre- 
quencies, and for oscillations with a given frequency, but 
for different excitations, can be used as a further observ- 
able for testing the system. The computation of ratios 
of amplitudes for different ft needs knowledge of F(f2). 
In principle, this can be measured from the width of the 
resonance. Ratios of amplitudes for fixed fl and differ- 
ent radial or angular dependence of the perturbation are 
independent of F(r2) and rrinim- In principle, amplitudes 
as in Fig. [5] can be measured both by dissipation (Imx) 
or by the amplitude of the cloud oscillations (Rex). 

We conclude that the variation of the external poten- 
tial with a certain frequency and radial and angular de- 
pendence allows for an excitation of individual modes. 
The oscillation can then be detected by imaging the den- 
sity profile or, hopefully, by new non-invasive methods. 
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E. Outlook on other systems 

We mentioned in the beginning of this section that 
we expect the features of the Bose gas to be generic for 
other systems of cold atoms. At nonzero temperature this 
requires the applicability of the two-fluid model which is 
related to a two-component order parameter. The latter 
can be written as (j){x) = |</)(a;)|e*^''^' at each space-time 
point. If 0{x) is varying slowly in space one can define 
the superfluid velocity as its gradient, Vs oc V0. Thus our 
description remains valid for systems which are described 
by a complex scalar field which acquires a non- vanishing 
expectation value. This includes in particular composite 
bosons like the atom-atom dimers on the BEC-side of the 
BEC-BCS crossover. For an imbalanced system it would 
then be interesting to include the conserved spin degrees 
of freedom of the remaining fermions as an additional 
macroscopic variable in the hydrodynamic equations. 

At zero temperature any system with hydrodynamic 
equations of motion might be considered. An example 
can be found in Ref. [52| where the collective oscilla- 
tion frequencies of a one-dimensional dipolar quantum 
gas have been calculated. 

The polytropic equation of state P cx /i""*"^ is of par- 
ticular interest for scale-invariant systems like ultracold 
Fermi gases at the unitary point where the scattering 
length a diverges. At zero temperature only the chemical 
potential remains as a dimensionful quantity, the pres- 
sure has to be proportional to /z^''"'"^'/^. Indeed, intro- 
ducing the Bertsch parameter ^ we can write the chem- 
ical potential at the unitary point as /i = S,Sf, where 
ep oc n'^/'^ is the Fermi energy density. Using dP = ndfj, 
we arrive at P{fJ.) oc ^('^+2)/2 as claimed. This time we 
know the equation of state exactly, but only for zero tem- 
perature. At the unitary point and for T = the oscil- 
lation frequencies obey the exact Eq. ([B]), with a = 3/2 
for a three-dimensional spherical trap. 

Around the unitary point and in the dilute, three- 
dimensional regime, the equation of state can be 
parametrized as [ssj 
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with an additional parameter Q implying a shift 

5uj% _ 256 C 1 

~ 5257r^/ci.(n(0))a 

of the breathing mode. 



(23) 



(24) 



III. COLLECTIVE MODES 
THEORY 



GENERAL 



In this section we derive the eigenvalue problem (jlip 
for collective modes at T > from Landau's two-fluid 
ideal hydrodynamics [13, S^. The heart of every hy- 
drodynamic description lies in the expansion of general 



constitutive equations in terms of derivatives of the con- 
tributing fields. The range of applicability of this de- 
scription has already been discussed in the introduction. 
The constitutive relations consist in conservation laws for 
particle number, momentum, energy and additional con- 
served quantities. The lowest order of this expansion is 
static equilibrium and thus thermodynamics. The next 
order is called ideal hydrodynamics and describes non- 
equilibrium processes without dissipation. The latter is 
included in higher orders of the expansion. 



A. The eigenvalue problem from two-fluid 
hydrodynamics 

Observations on quantum fluids suggest that there is 
a critical temperature Tc such that for < T < Tc the 
macroscopic hydrodynamic motion can be divided into 
two kind of flows, a superfluid (subscript s) and a nor- 
mal fluid (subscript n) one. Their velocity fields Vs{x,t) 
and Vn {x, t) , respectively, enter the momentum density 
g = UsVs + UnVn with coefficients satisfying n = -|- n„, 
where n(x, t) is the density of the gas or fluid under con- 
sideration. We might be lead to the conclusion that 
the substance is built out of two different kinds of flu- 
ids. However, the true nature of superfluidity consists 
in quantum physics and our description is only formally 
correct. The two flows can be distinguished by the follow- 
ing properties. The superfluid velocity field is a gradient 
field and hence irrotational. Furthermore, entropy is only 
carried by the normal fluid part. The full hydrodynamic 
equations under these constraints were set up by Lev. D. 
Landau in a famous and still very worth reading paper in 
1941 [l^]. We concentrate on ideal hydrodynamics where 
dissipation terms are neglected and entropy is conserved. 
In the homogeneous case, the underlying conservation 
laws read 



— + div(g) = 0, 



— + dfeH^fc 

ot 



0, 



ds , / ^ N 

— + div(st>„) = 0, 



dt 



/i 







(25) 
(26) 
(27) 
(28) 



with stress tensor n^fe = nnVn,iVn,k+nsVs,iVs,k+PSik, en- 
tropy density s, pressure P = —e+Ts-\-^n+nn{vn — VsY, 
chemical potential and internal energy density in the 
rest frame of the gas, e. These equations have to be 
supplemented by the equation of state -P(/i, T) and the 
normal fluid density nn{iJ.,T). For dissipative hydrody- 
namics one would replace Eq. ^7} by the conservation 
of energy and include an additional term in Eq. (|28l) . We 
use dimensionless units described in appendix [XI setting 
h — kg — m — 1. 

In order to understand what changes when applying an 
external held, we flrst derive the local density approxima- 
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tion from a thermodynamic point of view. For this pur- 
pose consider a gas contained in a potential V^xt (x) such 
that local equilibrium is reached at each point of space. If 
we take two neighboring small but yet macroscopic parts 
of the gas their energies Ei and E2 and particle numbers 
Ni and N2 will adjust in a manner maximizing the en- 
tropy S — Si + S2 under the constraints Ei + E2 = E 
and A^i -I- N2 — N, respectively. This implies tempera- 
ture and the full chemical potential to be constant inside 
the trap. However, the full chemical potential is given 
by the Gibbs free energy per particle and thus reads 
{P{x),T) + Vext{x), where {P,T) 
is the equilibrium chemical potential of the homogeneous 
system as a function of pressure P and temperature T. 
We conclude that a system where the chemical potential 
fj, in the homogeneous equilibrium functions is substi- 
tuted according to PhomifJ^T) PhomifJ. - Vext,T) etc. 
behaves like a system trapped in an external potential 
of large spatial extend. This rule is called local density 
approximation. 

An external potential can be included in Eqs. (|26p and 
by virtue of force terms, 

^ + dkll^k ^ -nd.Vext, (29) 
ot 



(30) 



The static solution denoted by a subscript gives the 
expected expressions ^o{x) + Vext{x) = p, and VPo(2?) + 
na{x)\/Vext{x) = of static equilibrium. Using the 
Gibbs-Duhem relation dP = sdT + ndn and thus 
VPq = soVTo + ngV/io every space point, we con- 
clude sqVTq — in equilibrium. We recognize the right 
thermodynamic behavior to emerge naturally. Pressure 
and densities of mass, entropy and energy are constant in 
time but space dependent, while temperature is constant 
all over the trap and the chemical potential follows the 
rule of local density approximation. The parameter p. 
can be adjusted in order to get a certain particle number 
N. The fluid velocities Vg and Vn vanish in equilibrium. 

The next step towards our eigenvalue problem consists 
in expanding the two- fluid equations in small fluctuations 
around their static equilibrium solution such that time- 
dependent quantities only appear linearly. We write 



fi{t, x) = ^o(^) + Suit, f), 
T{t,x) = T + ST{t,x), 

Vs.n(t,x) ^ SVs,n{t,x). 



(31) 
(32) 
(33) 



In the same way we linearize P = Pq + SP, n — uq + 
Sn,s ^ So + Ss. Eqs. ([13, (gZ]), ^ and §UI then 
become 



dtSn + div{Sg) — 0, 

dtSvs + V5^i = 0, 

dt5g + V5P + 5nWext = 0, 

dfSs + d\v(so5vn) = 



(34) 
(35) 
(36) 
(37) 



with 5g = Us^qSvs + n-„^Qdvn- Using the relation VPq = 
soVTo+rioV/io we get V5P = soV5T-|-^nV/io+rtoV(5/i = 
sqVST — SnWext + n^S/S^,, since VTq — and V/ig = 
—Wext- This can be used to cast Eq. into the form 



nn.ndt {6Vn - SVs) + S0V6T = 0. 



(38) 



In the following we will consider a spherically symmet- 
ric harmonic trapping potential Vextir) — y^q^'^ ■ The 
generalization of our formalism to other radial potentials 
Vexti'^) = Vext{r) is straightforward. In Eq. ([50)) we give 
a formulation of the eigenvalue problem which holds for 
arbitrary external potentials Vext (x) and thus also allows 
for a treatment of anisotropic trapping geometries. 

We assume the fluctuations of the physical quantities 
to be a harmonic oscillation in time with frequency w. 



6^{x, t) 



^Sfi{x), 



(39) 



and similar expressions for the other time-dependent 
functions. The remaining spatial deviations are again de- 
noted by (S/i etc. Since time has completely dropped out 
of our system of partial differential equations, this slight 
abuse of notation hopefully does not lead to confusion. 
We recognize our obtained set of equations 

iioSn — div{nsflSvs + «n,o<5wri), (40) 
iujSvs = V5fi, (41) 



So 



W6T, 



nn,o 

iuiSs — drv{sQSvn) 



(42) 
(43) 



to be an eigenvalue problem for a differential operator 
acting on (5/x, ST, 5vs and Svn — Svs- 



B. Zero temperature 

We derive the zero temperature eigenvalue problem in 
terms of an equation of state given in the form P(/^). We 
introduce the notation fioir) = /i(0) — Vextir) = ^J.{0) — 
ujQr'^ /2. Since we have set h = ks = m — 1, we can 
express all dimensionful quantities with respect to luq, see 
App. [X] In the following we set luq = 1 m all formulas, 
which means that we are measuring time in units of the 
inverse trapping frequency. Thus, /^o(^) = m(0) — 
Furthermore, we define 



PSir) ^ (^f ) (..(.•)). Pr(r) ^ (0 



Mr)), 



where /i(0) is the chemical potential in the center of the 
trap. We can write Sn{x) — Pq^ [r)5^{x) because dP = 
nd/i. There are no temperature fluctuations present at 
T = and thus the only degrees of freedom are J/x and 
5vs described by Eqs. (HUl) and (HT|) . The latter equations 
can be decoupled. Writing no(r) = Po{r) we arrive at 
the Stringari wave equation [16] 



Lo^P^^S^i + div(Po^V(5/i) 



0. 



(45) 
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So far, Sfj, depends on the d-dimensional spatial coordi- 
nate X. Since we are dealing with a spherically symmetric 
trap we can classify the possible solutions to the eigen- 
value problem via an ansatz 



(46) 



where fim are spherical harmonics. (For a definition of 
spherical harmonics in d < 3 dimensions we refer to App. 
[Cl Note that we have I = 0, 1 for d = 1.) When applying 
this ansatz to Eq. (H5|) we benefit from the facts that 
V/zq • y fim = and drPinoii')) — —fPai^)- In addition 
we use the relations 



e, • ^5^l{x) = ( -g'{r) + -g{r) ] h 



A5ii{x)=[-g"{r) 



r 

21 + d- 



(47) 
(48) 



where eV denotes the unit vector pointing in radial direc- 
tion. We arrive at 



g'{r) r'flra. 



= uj'^g{r) - rg'{r) - Igir) 



-g'{r)\. (49) 



When substituting z — the eigenvalue problem is given 
by 



Ag{z) = i^^g{z) 
for g{z) = g{r) and the differential operator 



(50) 



P^^(z) V dz^ ^ ' dz 



2z-^ + l 
oz 



(51) 

This equation has to be fulfilled on the interval z G [0, i?^] 
where R is the radius of the static cloud defined by no{r = 
R) = 0. We do not allow for fluctuations to change the 
cloud radius since this would be a second order small 
effect. Since we have set to 1, only appears in Eq. 
([50)) instead of w^/wq. 



C. Non- vanishing temperature 

For temperatures T > we have to implement tem- 
perature fluctuations and the full thermodynamic func- 
tions P(/i,T) and n„(/i,T). Extending the notation in- 
troduced in Eq. (|45l) we write 



(Mo(r),T). 



(52) 



In the same manner Pq, Pq'^^Pq, Pq'^, uq, and Uq are 
defined for T > 0, where 



The static solutions no and sq are connected to /xo(r) = 
- 7-2/2 and Tq = T via uq = P^, sq = Pq ■ Working 
again with independent variables Sfi and ST we get 

Sn = P^^'Sfi + PfST, 



Ss^P^^'Sh + Pq^'ST. 



TT 1 



(54) 
(55) 



From our analysis at zero temperature we know that due 
to the spherically symmetry the eigenmodes are most eas- 
ily obtained by eliminating the velocity fields. We find 



uP'Sn 



+ div(noV(5^ + sqV5T) = 0, (56) 
uj'^Ss + div(soV5M + noV<5T) = 0, (57) 



or in terms of P{fi,T), 
uj^P^'^Sfi + uj^P^'^ST + diY{P^V6fi + Pq^VJT) = 0, 



(58) 



oj^Po'^Sfi + i^j'^P^ ^' ST + div(P^ V(5/i + noVST) = 



(59) 



determining w. Note that no supplements the equation of 
state in the latter equation. We have chosen no instead 
of n„o in order to keep the notation simple. 

Eqs. ([58]) and (|59| can be written in the canonical form 
of an eigenvalue problem. Indeed, they are equivalent to 



pt^T pTT 
^0 ^0 



-div(Po^V.) 
-div(Po^V.) 



-div(Po^V- 
-div(noV-) 



ST 



Sfi 
ST 
(60) 



where we assumed Pq^ — P^^ . The first matrix appear- 
ing in this equation is the inverse of the Hessian matrix of 
P{fi,T). The latter is positive for thermodynamic rea- 
sons and thus the inverse always exists. The external 
trapping potential Vc^t{x) enters Eq. (|60l) through the 
application of the local density approximation indicated 
by the subscript 0. 

Analogously to the zero temperature case we clas- 
sify the solutions for an isotropic harmonic potential 
Vcvit{x) = /2 via an ansatz 



Sfj.{x) = gir^yji, 
ST{x) = Hr'^yf, 



(61) 
(62) 



where I — 0, 1, 2, ... and fim are the corresponding spheri- 
cal harmonics. The resulting equations simplify on going 



over to the new variable z 



We arrive at 



(53) 



A B 
C D 



ff(^) 
h{z) 



9{z) 
h{z) 



(63) 
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where the operators A, B, C, D are defined by 
A =«i(-) (4-^ + 2(2/ + rf)^) + (2^^ + /) (64) 
B ^h{z) Uz^ + 2(2/ + d)-^) + h{z) (2z-^^ + I 



C=ci{z){iz-^ + 2{2l + d)^_ 



dz^ 
d 



(65) 
(66) 



D =di(z) 4z— + 2(2/ + d)^ +d2{z) 2z— + l 



d 



with z-dependent coefficients 



d 



dz 



(67) 



ai 



pT pTT pM 



vt^T- pTT pT 



detPo 



detPo 



(68) 



pTT pA'T _ pf^T-fj. 
-'0 ^0 -'0 "-0 

dctPo 



, Ci 



J^O -Tq -Tq no 

detPo 



detPo 

(69) 

pMM^iM (ptJ.T\2 



detPo 



(70) 



with detPo = Po'^PP - (^o''^)^- Again, w = VI is a 
solution for constant g{z) and vanishing h{z), which is 
independent of the equation of state. 

An additional complication to the case of vanishing 
temperature arises from the fact that the normal part 
of the background density n„o is not restricted to a re- 
gion z < Zmax as it is the case for the superfluid density 
riso- In contrast it decreases exponentially for large z 
(c.f. Fig. HI). In praxis, we expect that the region z > E? 
for a sufficiently large value of is not affecting the 
oscillation frequencies since only an exponentially small 
number of particles is in that region. We restrict our 
numerical treatment therefore to a finite sphere z < R^. 

Given ai,...,d2 over the whole range z g [0, P^] the 
problem of determining the collective frequencies is as 
straightforward as it was for T — 0. Indeed, as we show 
in Sec. IIIIDI a simple and yet efficient method to solve 
for the eigenvalues will be given by replacing the opera- 
tors by matrices. Thus, if the more complicated coeffi- 
cient functions - (ffOl) are given to a sufficient degree 
of accuracy, the finite temperature collective modes are 
easily obtained. 



D. Numerical implementation 

Now that we have identified the ordinary differential 
operator(s) describing two-fluid hydrodynamic modes, 
we present a numerical method to determine the corre- 
sponding eigenfrequencies. The input for our approach 
consists in the equation of state P(/i, T) and normal fluid 



density nn{fJ,,T) on the interval G [^min, A*(0)]. Since 
we restricted ourselves to the case of spherically symmet- 
ric trapping the partial differential hydrodynamic equa- 
tions in d variables xi, . . . ,Xd are projected onto an or- 
dinary differential equation in z. This is a disadvantage 
of our method when compared to the fact that most ex- 
periments are performed in asymmetric traps. 

In order to introduce our numerical procedure we con- 
sider the case of T = since the operators for T > have 
an identical structure but different coefficient functions. 
We solve the eigenvalue problem (jH), ([5|) by discretiz- 
ing the interval z € [0,Zmax] via i — 0,1,..., M and 
Az — Zmax/M such that Zi — iAz. The function g{z) 
becomes a vector g — {g{zi))i=o^,,,^M and A is represented 
by a matrix (Aij) having components 



A, 



2z. 



S^- 



2A2 



-15--^^ 



X 4z, 



A^ 



2(2/ + <i) 



2Az 



(71) 



for i = 1, M — 1 and j = 0, M with Sij being the 
Kronecker delta. The cases i = and i = M require 
some care. We define 



oj ='ooj - pW(o) 2(2/ + d) 



Az 

/^A/j- 



(72) 



-^^r^2(2/ + d)^^^^^-^ (73) 

having replaced the second-order difference scheme by a 
first order one and set the second derivative g" at the 
boundary z = Zmax to zero. The eigenvalues of (Aij) can 
now be determined using standard procedures. 

Our discretization of A is closely related to the un- 
derlying physical problem. The functions g(z) on which 
A operates describe macroscopic hydrodynamic fluctua- 
tions of the chemical potential. Thus they are implic- 
itly assumed to be finite and free of a rich local struc- 
ture including jumps and kinks. This applies in partic- 
ular to the lowest lying modes with only a few radial 
nodes. Therefore we can safely assume g{z) to be suf- 
ficiently smooth such that the approximations g'{zi) = 
{g^+l~g^-l)/{2Az) + 0{Az^) and g"{z,) = {g,+i-2g,+ 
(7i_i)/(Az^)-|-0(Az^) hold in the open interval {O^Zmax)- 
However, at the points zq = and z^/ = z^ax we run into 
problems as z_i and zm+i are not defined. A first guess 
might be to assign certain boundary values to g. But 
investigating the breathing mode for a = 1 and d = 3 in 
Eq. (jD6|) . g[z) = 1 — 5z/6/i(0) for z < Zmax, otherwise, 
we recognize none of the sensible boundary conditions 
g{zmax) = or g'{zmax) = to be satisfied. This does 
not come as a surprise since our description necessarily 
breaks down in the outer regions of the cloud where the 
static density goes to zero. The crucial point is to pre- 
vent the solution from diverging at the boundaries. So we 
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(i) exact 


(i) num. 


(ii) exact 


(ii) num. 


(iii) exact 


(iii) num. 


3.46410 


3.46410 


1.50214 


1.50214 


1.00000 


1.00000 


8.00000 


7.99982 


2.35339 


2.35337 


2.16025 


2.16025 


12.4900 


12.4894 


3.13786 


3.13777 


3.10913 


3.10913 


16.6706 


16.9693 


3.89609 


3.89591 


4.00000 


4.00000 


21.4476 


21.4452 


4.64095 


4.64064 


4.86484 


4.86484 


25.9230 


25.9191 


5.37802 


5.37752 


5.71548 


5.71548 


30.3974 


30.3916 


6.11010 


6.10936 


6.55744 


6.55744 


34.8712 


34.8631 


6.83880 


6.83775 


7.39369 


7.39369 


39.3446 


39.3335 


7.56510 


7.56367 


8.22598 


8.22598 


43.8178 


43.8031 


8.28963 


8.28773 


9.05539 


9.05539 



TABLE I: Exact and numerical frequencies in units of ojo 
obtained for an equation of state P — ^""""^ for (i) a — 0.1, 
d = 1, I = Q, (ii) a = 3.9, d = 1, I = Q and (iii) a = 3, d = 3, 
1 = 1. All results correspond to a grid of M = 2000 points. 



where 1j\^{z) = jhy is the polylogarithm. Our nu- 

merical solution on a finite grid requires an effective ra- 
dius of the cloud although this is not a physical quantity. 
However, for too high distances from the center of the 
trap the gas will be so dilute that local equilibration can- 
not be ensured and therefore hydrodynamics does not 
provide the correct equations of motion for this part of 
the cloud. We therefore define an effective radius or min- 
imal chemical potential such that n(/Xmin)^o = 1' where 
is the oscillator length. At this value of the chemical 
potential the interparticle spacing will be of the order of 
the oscillator length and this constitutes the crossover re- 
gion where hydrodynamics breaks down. Since only very 
few atoms are in the outer regions of the cloud we expect 
the experimental hydrodynamic modes to be close to the 
modes which are calculated with this effective radius. 



demand it to be Taylor expandable even there, although 
the physical solution being zero for z > z^ax might have 
a discontinuity, and replace the second-order difference 
scheme in z_i and zm+i by a first-order one, i.e. for 
^0 = 

^i-^ = ^ + 0(Az). (74) 
2Az Az 

This can be solved for go and yields go ~ (gi + (7_i)/2. 
The latter mean-value property is true for any func- 
tion which can be linearized around z — provided 
a sufficiently small step size. In the same way, the 
approximation of g'{zmax) in Eq. (|73l) can be justi- 
fied. Setting g"{z„iax) = in the same equation implies 
dM+i — '^9m + 9m -1 = in a discretized version, which 
leads us to the same mean-value property. 

This simple approach is very efficient. As an example 
we show in Tab. U the first ten eigenfrequencies corre- 
sponding to a polytropic equation of state P{p) oc ^"^^ 
for different choices oi a, d and I compared to the exact 
solution given by Eq. ({6]). The underlying grid of points 
is of size M = 2000. We recognize the (more interesting) 
lower lying frequencies to converge faster. Of course, the 
accuracy can be improved by further increasing M. 

For a particular physical situation it may be favorable 
to use a non-uniform grid, i. e. the points Zi shall not 
be equidistant. Of course, the discretized operator (A^ ) 
can also be derived for a non-uniform grid. Such a dis- 
cretization has been applied to calculate the data points 
for Fig. [5] The use of a finite grid can be accompanied 
by a systematic error which overestimates or underesti- 
mates the eigenfrequencies. Nevertheless, we believe our 
method to be reliable because the exact results in Eqs. 
([6]) and (fTTl) are obtained to a sufficient accuracy and we 
found a quite fast convergence of the lowest eigenvalues 
from smaller to larger matrices. 

For a Bose gas at nonzero temperature the outer re- 
gions of the cloud are in the normal phase and the den- 
sity decays exponentially for large radii, n ~ Li3/2(e^^), 



IV. CONCLUSIONS AND OUTLOOK 

In this paper we suggest that collective oscillations may 
ultimately be used as precision observables for trapped 
ultracold atom gases. Within the hydrodynamic approx- 
imation the frequencies of these oscillations directly re- 
flect the thermodynamic equation of state. Their mea- 
surement can therefore shed light on the validity of non- 
perturbative methods for interacting many-body systems 
which are used for the computation of the equation of 
state. In parallel, further experimental advances may al- 
low for a simultaneous measurement of several frequen- 
cies, which could permit the determination of thermody- 
namic variables, in particular the temperature, with high 
precision. 

In the two-fluid hydrodynamic regime, the excitation 
and measurement of collective oscillations results in a 
set of discrete numbers which can be used to compare 
different theoretical methods. An experimental improve- 
ment of the precision of the measurement of collective 
modes would set stringent bounds on predictions from 
many-body calculations. Indeed, we have shown that 
the equation of state allows us to determine the oscilla- 
tion frequencies as eigenvalues of a differential operator. 
Thus, if broad collections of data on collective oscilla- 
tion frequencies would exist, this would necessarily rule 
out certain equations of state and theoretical approaches. 
In a certain sense this may be compared to restrictions 
of the composition and equation of state of the sun by 
helioseismology. 

The dependence of the collective modes on system pa- 
rameters as the scattering length is of particular interest. 
We have seen that interaction effects can lead to shifts 
of the zero temperature mean fleld results for the oscil- 
lation frequencies of a dilute Bose gas. For the three- 
dimensional case these shifts are rather small, however. 
Nevertheless, even there non-perturbative effects can be 
made visible. The appearance of an effective UV-cutoff 
related to the s-wave scattering approximation manifests 
itself in a dependence of the equation of state on the gas 
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parameter and this cutoff. Beyond lowest order perturba- 
tion theory we describe the effects of the LHY-correction 
for the oscillation spectrum. We also explore the range 
of larger scattering length, where perturbation theory is 
no longer valid, by means of an equation of state from 
Functional Renormalization Group (FRG). For gas pa- 
rameters n(0)a^ above 10~^ the difference between vari- 
ous methods for the relative frequency shift reaches the 
percent level and may be measurable. 

The analysis of the dilute Bose gas in two dimensions 
revealed further interesting many-body effects. Using 
again the equation of state from FRG there are sub- 
stantial deviations of the collective modes from mean 
field theory for large values of the dimensionless coupling 
constant. In these cases the difference of the breathing 
mode from its mean field value was found to be of order 
1% — 10% and showed an extremum for a specific value 
of the chemical potential in the center of the trap. Very 
interesting would be the investigation of these effects in 
a dimensional crossover from three to two dimensions. 

A main result of this paper concerns the temperature 
dependence of the oscillation frequencies. We suggest 
that this could ultimately be used as a precision ther- 
mometer. As an example we have computed the temper- 
ature dependent oscillation spectrum based on an equa- 
tion of state obtained from Lee Yang theory. The collec- 
tive oscillations of the condensate and the thermal cloud 
show a rich spectrum with non-trivial dependence on the 
temperature, including level-crossing and the disappear- 
ance of the superfluid modes at the critical temperature. 

The treatment of the oscillation spectrum at non- 
vanishing temperature has to account for the coexistence 
of two interacting fluids, one for the superfluid conden- 
sate and the other for the thermal cloud of uncondensed 
atoms. The thermal cloud disappears for T — ^ 0, while 
the condensate vanishes as the critical temperature is ap- 
proached, T — > Tc- Any computation of the spectrum has 
to reproduce these limiting cases correctly. In particular, 
we show how the zero temperature limit of the two-fluid 
hydrodynamic equations has a solution corresponding to 
a zero temperature formula for the condensate oscilla- 
tions given by Stringari. We also derived a new analytic 
result for frequencies at very low temperatures which are 
carried by the components of the thermal cloud. Our 
method of calculation covers the critical temperature as 
well as T > Tc, where the thermal cloud is the only fluid 
component. It will be interesting to study the effects 
of critical physics related to the second order superfluid 
phase transition on the oscillation frequencies. 

In order to measure not only the lowest collective os- 
cillation frequencies but rather parts of the whole spec- 
trum we suggest to use response techniques and excite 
frequencies individually. Once not only the lowest lying 
frequencies are accessible within an appropriate accuracy, 
but also higher modes, these will allow for precise verifi- 
cations of the consistency of theoretical models. 

The formal developments of this paper concern a 
method for the computation of oscillation frequencies 



in a trap, given the equation of state and normal fluid 
density for a homogeneous liquid. Starting from Lan- 
dau's ideal two-fluid hydrodynamics we have derived the 
general expression for a differential operator whose in- 
put is the equation of state given as pressure P(/x,T) 
and normal fluid density n„(^,T). The eigenvalues of 
this operator are the collective oscillation frequencies 
of the trapped gas. For simplicity we restricted our- 
selves to the case of an isotropic, d-dimensional, harmonic 
trapping. While the extension to an arbitrary spheri- 
cally symmetric trapping potential VextiT) is straightfor- 
ward, anisotropy would invalidate our symmetry assump- 
tions. Experiments are often performed in not completely 
spherically symmetric arrangements. We expect our find- 
ings to be qualitatively correct even for these cases. We 
propose a method for numerical implementation which 
is based on discretization. This can in principle also be 
done for an anisotropic trapping potential but with a 
higher computational effort. 

The restriction to the case of ideal hydrodynamics can 
be released. To first order, our equations would be mod- 
ified due to the five transport coefficients of superfluid 
hydrodynamics, including shear viscosity and heat con- 
ductivity. Here, immediately many interesting questions 
arise: How are these transport coefficients related to the 
damping rates of the modes, which appear as the phe- 
nomenological introduced widths in our treatment of lin- 
ear response? Can we use collective modes for measure- 
ments on the transport coefficients, e.g. rj/sl Is it pos- 
sible to derive, say, r]{fi,T) from a microscopic quantum 
field theory just like the equation of state T)l 

We conclude that collective modes have a high poten- 
tial to be used as precision observables, which can quan- 
tify different aspects of many-body physics. It seems 
worth to improve the techniques used for a precise mea- 
surement of oscillation frequencies. In particular, a non- 
invasive simultaneous measurement of several frequencies 
may allow for a precision estimate of the thermodynamic 
variables temperature and density or chemical potential. 

Appendix A: System of units 

In a relativistic system one usually sets h = ks = c = 1 
in order to express all quantities in terms of energies. 
However, for a non-relativistic system the velocity of light 
c does not appear in the equations and it is preferable to 
set h = ks = m = I, where m is a mass. (In our case 
this is the mass of the gas atoms.) All quantities can 
now be expressed in terms of a length scale of choice. 
The spherically symmetric trapping potential 

Vextix) = ^Ujy (Al) 

with = Xi + ■ ■ ■ + x-^ provides such a scale given by 
the oscillator length £o = ^/h/muo. If all dimensionful 

quantities are expressed in terms of the latter will no 
longer appear in the equations. In other words, we can 
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set £o = 1 in our formulas and understand all quanti- 
ties to be measured in units of Iq. Since h — m = 1 
this is equivalent to expressing everything in powers of 
ujQ. Thus, the eigenvalues of Ag{z) = uPg{z) have to be 
understood as dimensionless quantities lo/luq. 

We remark that due to m = 1 mass density equals 
particle number density, p — n, which simplifies our hy- 
drodynamic equations. 



Appendix B: Lee— Yang theory for Bose gas 

The equation of state for a weakly interacting Bose gas 
in three dimensions has been calculated by Lee and Yang 
in Refs. [soj . Starting from the microscopic dispersion 
relation ujp they compute the canonical partition function 
and derive the free energy density and other thermody- 
namic observables. The system is found to have two dis- 
tinct phases, a normal and a superfluid one. Besides this 
Lee and Yang also deduce the two-fluid hydrodynamics 
for the weakly interacting Bose gas. These equations in- 
clude that the mass m* of the collective excitations of 
the normal fluid component does not have to coincide 
with the mass m of the atoms. For the present work we 
neglect such details and only give a short summary of 
the equation of state which enters the Landau two-fluid 
hydrodynamics. 

Let density n and temperature T be given. We then 
have to determine the normal fluid density ri„(n, T) such 
that 



1 



where / = Jd^p/{2n f/^, £p = /2 and 



+ 2g{n - n„)) 



(Bl) 



(B2) 



is the microscopic dispersion. Since the latter depends on 
n„ itself, Eq. (|Bip has to be solved self-consistently. For 
zero temperature, where the normal fluid fraction van- 
ishes, we arrive at the standard Bogoliubov dispersion 
relation {sp{ep + 2gn)y^-^. We introduced a second func- 
tion j2{n,T) and thus we have to provide an additional 
equation. For a critical density 



n,{T) = C(3/2) 



3/2 



(B3) 



this equation reads Un — n in the case of n < ric, which 
corresponds to the normal phase. For n > ric we are in 
the superfluid phase and have 



1 



(B4) 



The (more complicated) latter case with two equations 
for Un and fl can be solved numerically by iteration. 



Once nn{n,T) and fl{n,T) are known the free energy 
density is found to be given by 



f{n,T) = |(n2 + nl) + f,n^ +T j log(l - e'^^'^--'^)). 



(B5) 



From /i — {df/dn)T we derive the chemical potential 



p. = il + g{n + n„). 



(B6) 



The pressure is then given by P{fi,T) = nfi — f, where 
n{fi,T) is found from inversion of Eq. (|B6p . Notice that 
the simple (exact) relation between n and fi for fixed n„ 
and /i allows us to use the grand variables /x and T right 
from the beginning. Then n„(/i,T) and fl{fi,T) have to 
be found self-consistently from Eqs. (|B4p - (IB6p using 
the dispersion relation 



y/£p{£p + 2(m - A) - 45n„), p. > fic 



(B7) 



instead of Eq. (|B2[) . The critical chemical potential is 
given by fic — 2gnc. 

The derivatives of the pressure can now be obtained in 
different ways. For the normal phase where the equations 
are particularly simple and involve only polylogarithmic 
functions one may use the canonical expressions derived 
from the free energy density and then insert n{fi, T) from 
Eq. (jB6p . For example from the pressure P(n, T) we can 
built {dP/dn)T = P^/P'^'^. Since P'^ = n is known, we 
get an expression for P^^ ^^l, T). Many expressions which 
are necessary for this process are given explicitly in the 
paper of Lee and Yang. For the superfluid phase the 
authors also give formulas for certain ranges of temper- 
ature. However, we cannot rely on these limiting cases 
since we want to cover the full equation of state for all 
values of /i and T . We found that the numerical imple- 
mentation of the integral expression of P(/x, T) derived 
from Eqs. (IB5[) and (|B6[) is suitable to calculate all other 
thermodynamic observables via finite difference formulas 
to sufficient accuracy. 



Appendix C: Spherical harmonics in one and two 
dimensions 

We define spherical harmonics in d dimensions for 
I — 0,1,2,... as complex valued functions fim on the 



unit sphere S"* ^ = {af G 



= 1} through 



the relation A{r fim) = 0, where A is the Laplacian, 



thonormality condition 



and demand the or- 



flmfl' 



(CI) 



m is a set of further indices. For d = 3 we have m G 
{—I, and fim = in traditional notation. In d = 2 
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dimensions, writing xi = rcos0, X2 = rsui(f>, we need to 
solve 

= A(rV,™(0)) -[^ + ~ + ^^) r'hM, 

(C2) 

(C3) 



which leads to 



Ji.±i = " 



/2tt 



for / = 0, 1, 2, ... We see that all values of I are allowed 
for d — 2. This is no longer true in the one-dimensional 
case as can already be guessed by physical reasoning. For 
d = 1 there are only monopole and dipole modes. Indeed, 
§° — {—1,1} is a discrete point set and one-dimensional 
spherical coordinates read x = ar with r = |a;| and a = 
x/r. The angular variable thus tells us whether we are 
considering the positive or the negative half line of the 
real axis. The expression g{x) = gi'f')fi(a) has to be 
understood as 



gir)M~l), x<0 
9ix) = { 5(0), X = 
g(r)/;(l), x>0. 



(C4) 



Since the harmonic functions in d — 1, i.e. the ones sat- 
isfying g" {x) — 0, are exactly the affine linear functions, 
the relation A(r'/;(o')) — can only be true for ^ = 0, 1. 
The normalization condition 



so 



is valid for 



Ilmfl'ni' - //*m(-l)//'m'(-l) + //m (1) /i'm' (1) 



^°^"^= 1/V2, .^1 ' 



(C5) 



cr = -1 
I/V2, a = 1 



(C6) 



/o and /i defined in this way have the correct parity 
(— 1)'. We do not give explicit constructions of spherical 
harmonics for higher dimensions than three. 



Appendix D: Exact results at zero temperature 

In this appendix we derive Eqs. (jB]) and (jl7p for col- 
lective modes in the zero temperature limit. These cal- 
culations provide an intuition on how solutions to our 
proposed eigenvalue problem should look like. Moreover, 
the exact eigenfrequencies allow to perform perturbation 
theory beyond mean field. In such a way Eq. (|9]) has 
been obtained in Ref. [s^. 

We start from the zero temperature eigenvalue problem 
Ag = uP'g with A from Eq. ([SJ . For a polytropic equation 



of state, P{\i) oc /x^+i, the ratio P^{z)IP^^{z) defined 
in Eq. (gl]) is found to be 



(q + 1)(a*(0)-z/2)° _ /i(0) - 
(a -t- l)a(/i(0) - z/2)°-i ~ < 



(Dl) 

where /i(0) is the chemical potential in the center of the 
trap. The eigenvalue problem then becomes 







L>?g{z) 



+ 2zg'(z) 
z/2 



a 



{Azg"{z) + 2{21- 



■d)g'{z)). (D2) 



This equation has to be true for all values of z inside the 
cloud. The radius of the cloud R is found by solving = 
no(r = i?) = (a + l)(^(0)-i?V2)", i.e. z < = 2^(0). 
Substituting z i— > z = z/2/i(0), the chemical potential 
drops out, 







uj'^giz) 
l~z 



2a 



2zg'(z) + lg{z) 
{A-zg'\z) + 2{2l + d)g'{z)). 



(D3) 



As is pointed out in Ref. [3] this is of the form of the 
hypergeometric difi'erential equation. However, for arbi- 
trary values oi a, d and I we would have to distinguish 
many special cases of its solution because the prefactors 
of the derivative terms in Eq. (|D3[) change. Therefore, 
we apply a more direct method and ask for solutions of 
the form g[z) 
relation 



- X]fe>o°^'=-^ ■ This yields the recursion 
u?) 4- 2k{a + fc + Z - 1 + rf/2) 



2(fc + / + rf/2)(fc + l) 



-a,. (D4) 



for the coefficients a^. Since Ag(z) — uP'g(z) is invariant 
under a rescaling of g{z) we can set — 1. A sufficient 
condition for the convergence of this series is termination 
and thus reduction to a polynomial of degree n = 0, 1, ... 
This will be the case for fc < n and 



2n 



1/2 

{a + n + l + d/2-l) + l\ ioo, (D5) 



where we restored the unit loq. Inserting this into Eq. 
(ID4p we get coefficients 



(fc - n){a + k + n + l + d/2-l) 
{k + l + d/2){k + l) 



-cik 



(D6) 



and Oo = 1. 

In the zero temperature limit the eigenvalue problem 
at nonzero temperature gets block diagonal because 
B and C vanish. While A approaches its zero temper- 
ature limit ([5|), the coefficient functions di and d2 of D 
approach 



m 

-pTT^ 



pTT 
^0 



(D7) 
(D8) 
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with no = (Pff)'^ /unfi, see Eq. ([53)) . For the dilute Bose 
gas in three dimensions these ratios can be computed 
analyticahy from the formulas given in App. |B] In the 
zero temperature limit nearly all values of the chemical 
potential obtained via ii{z) — /i(0) — Vext{z) satisfy the 
condition T <^ fJ-{z). This will of course not be true 
for /i = but this value can be neglected because it 
corresponds to the outermost points of the cloud. For 
T <^ fi the integrals for P and rin receive contributions 
only from the phonon part of the spectrum. We have 



(D9) 
(DIO) 



see for example 0. These quantities are powers of T 
and will be extremely small. Nevertheless, the ratios 
(|D7p and (jPSp arc non-vanishing because the powers of 
T cancel each other. We find 



(Dll) 
(D12) 



and the differential equation D{T — 0)h{z) — u)'^h{z) for 
thermal fluctuations becomes 

DiT = 0)h{z) = - -(2zh'(z) + Ihiz)) 
6 

- i^^^h-lH{Azh"{z) + 2{2l + i)h'{z)) 



io^h{z). 



(D13) 



We already inserted d = 3. When we multiply both sides 
with —6 we find the same structure as in Eq. (jD2l) but 
with the substitutions oj^ —60;^ and a n- —1/2. Ap- 
plying these modifications to the square of both sides of 
Eq. (lD5l) . 



2n 
a 



a + n + l + l: ) +1 



we find the eigenvalues of Eq. (|D13|) to be 

' An{n + I) ~ iV^^ 
w„,i = 1 I Wq. 



(D14) 



(D15) 



This proves Eq. (|17l) . The eigenfunctions read h{z) — 
Sfc=o'^fe^'^ with z — z/2/i(0), ao = 1 and 



__ k{k + 1) - n{n + I) 
= {k + l + 3/2){k+lf'- 



Appendix E: Response to an external driving force 



In addition to the trapping potential Vext{x) — ^uj^r^ 
we want to apply a small perturbation SVext {x, t) such 
that the oscillating system experiences a driving force. 



We will show in this appendix how the equations of a 
driven oscillator arise and which perturbations can be 
used to excite the eigenfrequencies and thus measure the 
spectrum of collective modes. In order to keep the deriva- 
tion short we only consider the zero temperature case. It 
is straightforward to generalize the following treatment 
to non-vanishing temperatures. 

When we include the perturbing potential the lin- 
earized equation ((35)) for the superfluid velocity becomes 
dtSvs + W{Sfi + SVext) = 0. Hence the Stringari wave 
equation ([i5|) is now replaced by 

{d^ + rdt)PrSt^-div{P^W6fi)^div{Pi;WSVext), (E1) 



where we also include damping like in Eq. (|18p . The 
left-hand side of this equation represents the unper- 
turbed oscillator and the right-hand side corresponds to 
the external driving force. Applying a Fourier transfor- 
mation Sii{x,t) = / ^e~*^*(5/i(af, £7) and analogous for 
SVext {x, t) we obtain 

{-n^ - ivvL + E)5fi{x, n) = ~ESVext{x, n), (E2) 

where we defined the operator E — — (P(^'')~^div(P^V-). 
Note that for 4>nimix) = gni{z)r^ fira we have E4)nim,{x) = 
r'- fimAgni{z) = w'^i4>nim{x) with A from Eq. ([H]) and 
corresponding eigenfunctions gni{z) of A. Thus, A and 
E share the same eigenvalues. 

We introduce the scalar product (/, 5) = 
J^f{x)*P^'^{r)g{x) where the integral extends over 
the cloud. The operator E is then self-adjoint and its 
eigenfunctions (/>„im form a complete and orthogonal 
system. For f = g = we get (/, g) = J SnS^, where 
the chemical potential is multiplied by its thermody- 
namic conjugate, the density. We assume to be 
normalized to 1 in the following. 

For a perturbing potential dVext (x, t) , the response 
function x is defined as 

6nix,t)^ f dT f xix,y,T)6Vextiy,t-T). (E3) 

Jo Jy 

Introducing x{x,y,^) = dte'^^'^x{x,y,t) this can be 
written as 

s^i{x,n)^ f x{x,y,n)dVext{y,n). (E4) 

The spatial dependence on the coordinate x can be elim- 
inated by expanding in the eigenfunctions (f>nim of E, 
which yields 



(D16) dflrU^M^ E 

Xnlm,n'l'm' {rt){SVext)n'l'm'{^)- (E5) 



n'l'm' 



From the equation of motion (jE2l) we can immediately 
read off the response function in this particular basis. It 
is given by 



\nlm ^n' I' m' iS^^ 



-Snn'Su'Smm'- (E6) 
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Eq. (jE6|) shows that the mode i 



can be excited 



by applying a periodic perturbation of exactly this spatial 
form. 

For practical purposes we are more interested in the 
response of the trapped gas to an arbitrary perturbation 
such that in general several modes will be excited. Let 
the shape of SVext be given by 



(E7) 



The variation of the external potential leads to a per- 
turbation jg5n{x)5Vext{x,t) of the Hamiltonian of the 
system. The rate of average energy dissipated is there- 
fore given by 



"dF 



— = ( / 5n{x)dt5Vext{x,t) 



(E8) 



where the averaging is over a small interval of time. Using 
5n = Pq^S^ this becomes 



— Im 
2 



H*{x)P^^{r)x{x,v,^l)H{y) 



(E9) 

For an oscillating system without spatial coordinate x, 
Eq. dUl would yield E = (17/2)|iJopInix(rj) for con- 
stant H{x) = iJg. If we consider the trapped system as a 
whole, we can use this analogy to define the response of 
the cloud to the perturbation H{x). It has a clear physi- 
cal meaning related to the dissipated energy. Expanding 
Hix) = J2nhn Hrdra(t)nim{x) wc arrive at 



dt 

nlm 

Therefore, we have 



n/m I 



(ElO) 



-1^2 - irn 



nl 



with response coefficients 
1 



— \Hnlm\ LUnl 



(E12) 



We simplify the situation even further by considering 
a three-dimensional mean-field Bose gas. As we have 
seen in Eq. this corresponds to an equation of state 
P(/i) = /i^/87ra. The eigenfunctions were found in App. 



|D]to be given by (j)nim{x) = gni{z)r^ fim with gni{z) = 
Z]fc=o"fc^'^' z = {t/RY and defined in Eq. (jD6p . 
The squared radius in this case is = 2/x(0) with the 
chemical potential in the center of the trap /i(0). The 
lowest isotropic modes {I = 0) are given by 



510 (^) 

93o{z) 
9io{z) 



14 21 . 

= 1 z-i 

3 5 

1 - 9z + —z^ - 



429 



44 



286 



572 



2431 



5 - 

We further assume H entering Eq. (jE7p to be of the 
form 



H{x) = Fir^fLM 



(E14) 



with monomial F{r) — r'^ and q > 0. We have 
(wnim)"^ SiL^mM- We coucludc that a perturbing 
potential of the form (|E14[) has to have the right angular 
behavior (oc /;,„) in order to excite a mode Sfj, <x fim- 
This is not restricted to the case of the Bose gas but 
holds in general. For L = Af = we find 



1 



miooiq) 
1 

m2oo{q) 
1 

msooiq) 
1 



2iq + 3)iq + 5)J 

165((7 - 2)^714 y 
8{q + 3){q + 5)(q + 7) ) ' 

525(g-4)(g-2)gV27 V 
16{q + 3){q + 5){q + 7){q + 9)) ' 

5985(g-6)(g-4)(g-2)gV44 V 

(Ell) m4oo(9) " [l28{q + 3)iq + 5)iq + 7)iq + 9)iq+n) ) 

(E15) 

Here, we omitted an overall prefactor of R'^^'^/Sttu. The 
response coefficients depend strongly on q, i.e. the mono- 
mial character of SVext ■ In particular for low values of q 
there will not be a response from the higher lying modes. 
For g = 0, i.e. 6Vext{t,x) = cos{nt), none of the collec- 
tive modes will be excited, for SVext = r^cos{^lt) only the 
breathing mode, and so on. We considered here isotropic 
oscillations. An analog behavior can be found for I > 0. 
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